#!/bin/bash

while getopts  ":s:p:g" opts
do
	case  $opts  in
	s) sample_name=$OPTARG;;
	p) out_prefix=$OPTARG;;
	g) summ_gatk=T;;
	esac
done
shift $(($OPTIND - 1))

if [ -z "$1" ]; then
	echo `basename $0` ' [-s sample_name] [-g] <input.bam> <interval.bed>'
	exit 65
fi


. $var


bam0=$1
# interval=$2

interval0=$out_prefix.bed
if test ! -e "$out_prefix.bed"; then 
cut -f1-3 $2 > $interval0
fi

if test "$summ_gatk" = "T"; then
summa.sh -p$out_prefix $1 $2 &
fi

depth=$out_prefix.depth.txt
depth0=depth.txt
sum=$out_prefix.sum.txt

m_path=${tools_path}/matlab


:> $sum
echo -e "Sample\t\t$sample_name" >> $sum

samtools flagstat $bam0 > $out_prefix.stat.txt

total_reads=`cat $out_prefix.stat.txt |  grep total | cut -d "+" -f 1`
paired=`cat $out_prefix.stat.txt | grep "paired in" |cut -d "+" -f 1`
paired_p=`echo "scale=2; $paired/$total_reads*100" |bc -l`
proper_paired=`cat $out_prefix.stat.txt | grep "properly paired" |cut -d "+" -f 1`
proper_paired_p=`echo "scale=2; $proper_paired/$total_reads*100" |bc -l`
mapped=`cat $out_prefix.stat.txt | grep -v mate|grep mapped |cut -d "+" -f 1`
mapped_p=`echo "scale=2; $mapped/$total_reads*100" |bc -l`
paired_mapped=`cat $out_prefix.stat.txt |grep 'mate mapped'|grep -v 'mate mapped to a different chr'|cut -d "+" -f 1`
paired_mapped_p=`echo "scale=2; $paired_mapped/$total_reads*100" |bc -l`
echo -e "total reads\t\t$total_reads" >> $sum
echo -e "paired reads\t$paired_p%\t$paired" >> $sum
echo -e "mapped reads\t$mapped_p%\t$mapped" >> $sum
echo -e "paired and mapped reads\t$paired_mapped_p%\t$paired_mapped" >> $sum
echo -e "properly paired reads\t$proper_paired_p%\t$proper_paired" >> $sum


# $bedtools genomecov -ibam 1.target.bam -bga -split -g $ref_genome > 1.bedGraph
echo;echo bedtools coverage
bedtools coverage -d -split -hist -abam $bam0  -b $interval0 > $depth
# $bedtools coverage -abam $bam0  -b $interval > $gene_depth

cut -f5 $depth > dpt
total=`sed -n '$=' dpt`
seq  $total > dpt.id
paste dpt.id dpt > $out_prefix.dpt
echo depth.py
depth.py -t$total -p$out_prefix.dpt $out_prefix.dpt 

if [ `pwd` != $m_path ]; then
	cp $m_path/summ.m $m_path/plot_depth.m .
fi

# awk '{print $6}' $depth > $depth0
awk '{print $5}' $depth > $depth0
# awk '{print $5}' $gene_depth > $gene_depth0

echo summ.m need the file $depth0 and output the file 1.summ.txt
matlab -nojvm -r "summ; quit;" -logfile .m.log


if [ `pwd` != $m_path ]; then
	rm summ.m plot_depth.m
fi

cat $sum 1.summ.txt > $out_prefix.summary.txt \
&& tab2xls.pl $out_prefix.summary.txt $out_prefix.align.cov.summary.xls

wait

rm $sum 1.summ.txt $out_prefix.stat.txt depth.txt dpt dpt.id $out_prefix.dpt $out_prefix.depth.txt

. $cmd_done


